{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "private_outputs": true,
   "provenance": []
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/mikexcohen/Statistics_book/blob/main/stats_ch12_correlation.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Modern statistics: Intuition, Math, Python, R\n",
    "## Mike X Cohen (sincxpress.com)\n",
    "### https://www.amazon.com/dp/B0CQRGWGLY\n",
    "#### Code for chapter 12 (correlation)\n",
    "\n",
    "---\n",
    "\n",
    "# About this code file:\n",
    "\n",
    "### This notebook will reproduce most of the figures in this chapter (some figures were made in Inkscape), and illustrate the statistical concepts explained in the text. The point of providing the code is not just for you to recreate the figures, but for you to modify, adapt, explore, and experiment with the code.\n",
    "\n",
    "### Solutions to all exercises are at the bottom of the notebook.\n",
    "\n",
    "#### This code was written in google-colab. The notebook may require some modifications if you use a different IDE."
   ],
   "metadata": {
    "id": "mzfXc9E3Xq7k"
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dFT1qwVEyxTW"
   },
   "outputs": [],
   "source": [
    "# import libraries and define global settings\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# import module from scipy (for cosine similarity)\n",
    "from scipy import spatial\n",
    "\n",
    "# pandas and seaborn for the exercises\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "# define global figure properties used for publication\n",
    "import matplotlib_inline.backend_inline\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats('svg') # display figures in vector format\n",
    "plt.rcParams.update({'font.size':14,             # font size\n",
    "                     'savefig.dpi':300,          # output resolution\n",
    "                     'axes.titlelocation':'left',# title location\n",
    "                     'axes.spines.right':False,  # remove axis bounding box\n",
    "                     'axes.spines.top':False,    # remove axis bounding box\n",
    "                     })"
   ]
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "AQV2yRQpbRco"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 12.1: Example of scatter plot showing correlated data"
   ],
   "metadata": {
    "id": "ChAn8kbBTIQr"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Number of sibling pairs\n",
    "num_pairs = 40\n",
    "\n",
    "# Simulate ratings for brothers\n",
    "brothers_ratings = np.random.randint(1,11,num_pairs)\n",
    "\n",
    "# Simulate correlated ratings for sisters based on brothers' ratings\n",
    "noise = np.random.normal(0,2,num_pairs)  # some random noise\n",
    "sisters_ratings = brothers_ratings + noise  # sister's ratings are brother's ratings plus some noise\n",
    "\n",
    "# Ensure ratings are within bounds 1 and 10\n",
    "sisters_ratings = np.clip(np.round(sisters_ratings),1,10)\n",
    "\n",
    "# correlation\n",
    "r = stats.pearsonr(brothers_ratings,sisters_ratings)[0]\n",
    "\n",
    "# Create a scatter plot of the data\n",
    "plt.figure(figsize=(6,5))\n",
    "plt.plot(brothers_ratings, sisters_ratings,'ks',markersize=10,markerfacecolor=(.7,.7,.7))\n",
    "plt.title(f'Death metal preferences in siblings (r={r:.2f})',loc='center')\n",
    "plt.xlabel(\"Brothers' ratings\")\n",
    "plt.ylabel(\"Sisters' ratings\")\n",
    "plt.xlim(0,11)\n",
    "plt.ylim(0,11)\n",
    "plt.grid(True)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_death.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "8jzfGd9ITIT3"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "PvFC7_5fTIWu"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 12.2: Different correlation coefficients"
   ],
   "metadata": {
    "id": "L4GyyIv3bRZ1"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# correlation values\n",
    "rs = [ 1,.7,.2,0,0,0,-.2,-.7,-1 ]\n",
    "\n",
    "# sample size\n",
    "N = 188\n",
    "\n",
    "\n",
    "# start the plotting!\n",
    "_,axs = plt.subplots(3,3,figsize=(8,8.5))\n",
    "\n",
    "for r,ax,i in zip(rs,axs.flatten(),range(9)):\n",
    "\n",
    "  # generate data\n",
    "  x = np.random.randn(N)\n",
    "  y = x*r + np.random.randn(N)*np.sqrt(1-r**2)\n",
    "\n",
    "  # exceptions for r=0\n",
    "  if i==3:\n",
    "    x = np.cos(np.linspace(0,2*np.pi-2*np.pi/N,N))\n",
    "    y = np.sin(np.linspace(0,2*np.pi-2*np.pi/N,N))\n",
    "  elif i==4:\n",
    "    x = np.linspace(-2,2,N)\n",
    "    y = x**2\n",
    "  elif i==5:\n",
    "    x = np.linspace(-2,2,N//2)\n",
    "    y = np.concatenate((x,-x),0)\n",
    "    x = np.concatenate((x,x),0)\n",
    "\n",
    "  # empirical correlation\n",
    "  rho = np.corrcoef(x,y)[0,1]\n",
    "\n",
    "  # plot\n",
    "  ax.plot(x,y,'ko',markersize=10,markerfacecolor=(.7,.7,.7),alpha=.3)\n",
    "  ax.set(xticks=[],yticks=[])\n",
    "  ax.set_title(f'r = {rho:.2f}',loc='center')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_variousRs.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "Epiwvc8JD7-X"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "T4oErk_Drtw8"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 12.3: Same correlation, different slopes"
   ],
   "metadata": {
    "id": "Sfnil-gOWbMo"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "N = 100\n",
    "\n",
    "# Dataset 1\n",
    "x1 = np.random.normal(100,10,N)\n",
    "y1 = .3*x1 + np.random.randn(N)*3\n",
    "slope1,intercept1,r1,_,_ = stats.linregress(x1,y1)\n",
    "\n",
    "# Dataset 2\n",
    "x2 = np.random.normal(10,1,N) + np.mean(x1)\n",
    "y2 = 3*x2 + np.random.randn(N)*3\n",
    "slope2,intercept2,r2,_,_ = stats.linregress(x2,y2)\n",
    "\n",
    "# x-axis limits\n",
    "xmin,xmax = np.min(x1)-5,np.max(x1)+5\n",
    "\n",
    "# Plot datasets and their regression lines\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "axs[0].plot(x1,y1,'ko',markersize=10,markerfacecolor=(.3,.3,.3),alpha=.3)\n",
    "axs[0].plot(x1,intercept1 + slope1*x1,'k')\n",
    "axs[0].set_title(fr'$\\bf{{A}}$) Slope={slope1:.2f}, r={r1:.2f}')\n",
    "axs[0].set(xlim=[xmin,xmax])\n",
    "\n",
    "axs[1].plot(x2,y2,'ks',markersize=10,markerfacecolor=(.3,.3,.3),alpha=.3)\n",
    "axs[1].plot(x2, intercept2 + slope2*x2,'k')\n",
    "axs[1].set_title(fr'$\\bf{{B}}$) Slope={slope2:.2f}, r={r2:.2f}')\n",
    "axs[1].set(xlim=[xmin,xmax])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_fitlineR.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "v7olQt-JWdCk"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "Rzgkt0uIMvzt"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# numpy vs scipy"
   ],
   "metadata": {
    "id": "NGDqeDOAMvxB"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# vector input\n",
    "x = np.random.randn(50)\n",
    "y = x + np.random.randn(len(x))/2\n",
    "\n",
    "usingScipy = stats.pearsonr(x,y)\n",
    "usingNumpy = np.corrcoef(x,y)\n",
    "\n",
    "print('scipy.stats.pearsonr:')\n",
    "print(usingScipy)\n",
    "\n",
    "print(' ')\n",
    "print('numpy.corrcoef:')\n",
    "print(usingNumpy)"
   ],
   "metadata": {
    "id": "ZRYePmtsMvul"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# matrix input (features by observations)\n",
    "X = np.vstack((x[None,:],y[None,:]))\n",
    "\n",
    "# usingScipy = stats.pearsonr(X) ## gives an error!\n",
    "usingNumpy = np.corrcoef(X)\n",
    "\n",
    "print('numpy.corrcoef:')\n",
    "print(usingNumpy)"
   ],
   "metadata": {
    "id": "l9GGmfI1Mvr4"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# to get a matrix of R and p values:\n",
    "\n",
    "R = np.zeros((2,2))\n",
    "P = np.zeros((2,2))\n",
    "\n",
    "# Calculate Pearson correlation for each pair of variables\n",
    "for i in range(2):\n",
    "  for j in range(2):\n",
    "    R[i,j], P[i,j] = stats.pearsonr(X[i,:],X[j,:])\n",
    "\n",
    "\n",
    "print('Correlation matrix:')\n",
    "print(R)\n",
    "\n",
    "print('')\n",
    "print('Associated p-values:')\n",
    "print(P)\n",
    "\n",
    "# Note about this code cell: you don't actually need to compute all matrix elements;\n",
    "# you can compute only the upper-triangle and then copy the results to the lower triangle."
   ],
   "metadata": {
    "id": "XpEgS-ktMvma"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "CT_1mPuNdxhb"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Creating correlated data"
   ],
   "metadata": {
    "id": "QzfiCPlFdxeQ"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Method 1 (quick&dirty but effective)\n",
    "x = np.random.randn(40)\n",
    "y = x + np.random.randn(len(x))\n",
    "r = stats.pearsonr(x,y).statistic\n",
    "\n",
    "plt.plot(x,y,'ks')\n",
    "plt.title(f'r = {r:.3f}',loc='center')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "k8OCxiKDd0M2"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# method 2 (more control over the correlation)\n",
    "r = .4\n",
    "x = np.random.randn(50)\n",
    "y = np.random.randn(len(x))\n",
    "y = x*r + y*np.sqrt(1-r**2)\n",
    "\n",
    "rr = stats.pearsonr(x,y).statistic\n",
    "\n",
    "plt.plot(x,y,'ks')\n",
    "plt.title(f'r = {rr:.3f}',loc='center')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "qq-_EML2dxbK"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# method 3: multivariate\n",
    "C = np.array([ [1,.4],[.4,1] ])\n",
    "means = np.zeros(2)\n",
    "\n",
    "X = np.random.multivariate_normal(means,C,50)\n",
    "r = np.corrcoef(X.T)\n",
    "\n",
    "plt.plot(X[:,0],X[:,1],'ks')\n",
    "plt.title(f'r = {r[0,1]:.3f}',loc='center')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "ISQCoMOHgW_D"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "rx56Vo8Uoq5t"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 12.6:  Anscobe's quartet"
   ],
   "metadata": {
    "id": "lSKQ1QHOK1P7"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "\n",
    "anscombe = np.array([\n",
    "     # series 1     series 2      series 3       series 4\n",
    "    [10,  8.04,    10,  9.14,    10,  7.46,      8,  6.58, ],\n",
    "    [ 8,  6.95,     8,  8.14,     8,  6.77,      8,  5.76, ],\n",
    "    [13,  7.58,    13,  8.76,    13, 12.74,      8,  7.71, ],\n",
    "    [ 9,  8.81,     9,  8.77,     9,  7.11,      8,  8.84, ],\n",
    "    [11,  8.33,    11,  9.26,    11,  7.81,      8,  8.47, ],\n",
    "    [14,  9.96,    14,  8.10,    14,  8.84,      8,  7.04, ],\n",
    "    [ 6,  7.24,     6,  6.13,     6,  6.08,      8,  5.25, ],\n",
    "    [ 4,  4.26,     4,  3.10,     4,  5.39,      8,  5.56, ],\n",
    "    [12, 10.84,    12,  9.13,    12,  8.15,      8,  7.91, ],\n",
    "    [ 7,  4.82,     7,  7.26,     7,  6.42,      8,  6.89, ],\n",
    "    [ 5,  5.68,     5,  4.74,     5,  5.73,     19, 12.50, ]\n",
    "    ])\n",
    "\n",
    "\n",
    "# plot data and correlations\n",
    "fig,ax = plt.subplots(2,2,figsize=(8,5))\n",
    "ax = ax.ravel()\n",
    "\n",
    "for i in range(4):\n",
    "\n",
    "  # plot the points\n",
    "  ax[i].plot(anscombe[:,i*2],anscombe[:,i*2+1],'ko',markersize=10,markerfacecolor=(.7,.7,.7))\n",
    "\n",
    "  # compute the corrs\n",
    "  corr_p = stats.pearsonr(anscombe[:,i*2],anscombe[:,i*2+1])[0]\n",
    "  corr_s = stats.spearmanr(anscombe[:,i*2],anscombe[:,i*2+1])[0]\n",
    "\n",
    "  # update the axis\n",
    "  ax[i].set(xticks=[],yticks=[])\n",
    "  ax[i].set_title(f'$r_p = {corr_p:.2f},  r_s = {corr_s:.2f}$',loc='center')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_anscobe.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "KgoVX3TbFiVi"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "qRXUYAtjbRV4"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Toy covariance example"
   ],
   "metadata": {
    "id": "gSKBJ-orbRVH"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# raw scores\n",
    "h = np.array([74,63,58,70])\n",
    "s = np.array([ 4, 7, 2, 9])\n",
    "N = len(h)\n",
    "\n",
    "# demeaned\n",
    "hd = h-np.mean(h)\n",
    "sd = s-np.mean(s)\n",
    "\n",
    "cov = np.sum(hd*sd) / (N-1)\n",
    "cov, hd*sd"
   ],
   "metadata": {
    "id": "Re0GDWeQGw5i"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "kInGUrv6bRM0"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 12.7: Kandall tau"
   ],
   "metadata": {
    "id": "kJrytL6QbRKF"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# The data\n",
    "bro = np.array([ 1,2,3,4,5 ])\n",
    "sis = np.array([ 2,1,4,5,3 ])\n",
    "\n",
    "# the correlation\n",
    "k = stats.kendalltau(bro,sis)\n",
    "\n",
    "# band names (lol)\n",
    "bands = [ 'Unicorn Apocalypse',\n",
    "          \"Satan's Fluffy Bunnies\",\n",
    "          'Demonic Dishwashers',\n",
    "          'Vampiric Vegetarians',\n",
    "          'Zombie Zucchini Zephyr' ]\n",
    "\n",
    "# the plot\n",
    "_,ax = plt.subplots(1,figsize=(6,5))\n",
    "ax.plot(bro,sis,'ks',markersize=30,markerfacecolor=(.9,.9,.9))\n",
    "ax.set(xlabel=\"Brother's ratings\",ylabel=\"Sister's ratings\",\n",
    "       yticks=range(1,6),xticks=range(1,6),\n",
    "       xlim=[.5,5.5],ylim=[.5,5.5])\n",
    "ax.grid()\n",
    "ax.set_title(fr'$\\tau$ = {k[0]:.2f}',loc='center')\n",
    "\n",
    "\n",
    "\n",
    "# band names\n",
    "for i in range(len(bro)):\n",
    "\n",
    "  # band names in plot markers\n",
    "  ax.text(bro[i],sis[i],bands[i][0],size=20,\n",
    "           weight='bold',va='center',ha='center')\n",
    "\n",
    "  # raw data\n",
    "  ax.text(9,i+1,bands[i],ha='right',va='center')\n",
    "  ax.text(9.5,i+1,bro[i],ha='center',va='center')\n",
    "  ax.text(10,i+1,sis[i],ha='center',va='center')\n",
    "\n",
    "# column labels\n",
    "ax.text(9,5.5,'Band name',ha='right',va='center',weight='bold')\n",
    "ax.text(9.5,5.5,'Bro',ha='center',va='center',weight='bold')\n",
    "ax.text(10,5.5,'Sis',ha='center',va='center',weight='bold')\n",
    "\n",
    "\n",
    "# plt.tight_layout()\n",
    "plt.savefig('cor_kendall.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "aCHEB4HNpDr-"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "9DpRMismbREk"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 12.8: Statistical significance of r based on n"
   ],
   "metadata": {
    "id": "CLFEBdKlbRBd"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# simulation ranges\n",
    "rs = np.linspace(.1,.8,53) # r values\n",
    "ns = np.arange(10,511,step=25) # sample sizes\n",
    "\n",
    "# compute the matrix of t-values\n",
    "num = rs[:,None]*np.sqrt(ns-2)\n",
    "den = 1-rs[:,None]**2\n",
    "tmat = num/den\n",
    "\n",
    "\n",
    "## show the matrix!\n",
    "fig,ax = plt.subplots(1,figsize=(4,6))\n",
    "\n",
    "cax = ax.imshow(tmat,vmin=2,vmax=10,aspect='auto',cmap='gray',\n",
    "           extent=[ns[0],ns[-1],rs[0],rs[-1]],origin='lower')\n",
    "ax.set(xlabel='Sample size',ylabel='Correlation strength')\n",
    "\n",
    "# and make it look a bit nicer\n",
    "fig.colorbar(cax,orientation='horizontal',pad=.02,ax=ax,location='top')\n",
    "ax.spines[['right','top']].set_visible(True)\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_tvals.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "SN2YeAWjbQ8m"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "rrtwur2mbhOd"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 12.9: Fisher-z transform on uniform data"
   ],
   "metadata": {
    "id": "LbXP4WCYS8YS"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# \"correlation\" data and its Fisher transform\n",
    "r = np.random.uniform(-1,1,size=1000)\n",
    "fish_r = np.arctanh(r)\n",
    "\n",
    "\n",
    "\n",
    "# now for plotting\n",
    "fig,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "# histograms\n",
    "axs[0].hist(fish_r, bins='fd', color=(.3,.3,.3),alpha=1, label='F(r)')\n",
    "axs[0].hist(r,bins='fd', color='k', alpha=.5, label='r')\n",
    "axs[0].legend(loc='upper right')\n",
    "axs[0].set_xlabel('Value')\n",
    "axs[0].set_ylabel('Frequency')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Histograms of r and F(r)')\n",
    "\n",
    "# scatter plot to visualize the transform\n",
    "axs[1].plot(r,fish_r,'ko',markersize=10,alpha=.5,markerfacecolor=(.8,.8,.8))\n",
    "axs[1].set_xlabel('r')\n",
    "axs[1].set_ylabel('F(r)')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Scatter plot of r vs. F(r)')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_fisherFull.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "xAG_LglpS_sO"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "wY8BZ4r51ngb"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 12.10: Fisher-z transform on numbers close to zero"
   ],
   "metadata": {
    "id": "yg6HL6Y71nuW"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Same as above but with simulated H0 coefficients\n",
    "\n",
    "# random data\n",
    "X = np.random.randn(100,45)\n",
    "\n",
    "# compute the correlation matrix and extract the unique r values\n",
    "R = np.corrcoef(X.T)\n",
    "utri = np.triu(R,k=1)\n",
    "r = utri[utri!=0]\n",
    "\n",
    "# fisher transform\n",
    "fish_r = np.arctanh(r)\n",
    "\n",
    "\n",
    "# now for plotting\n",
    "fig,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "# histograms\n",
    "axs[0].hist(fish_r, bins='fd', color=(.3,.3,.3),alpha=1, label='F(r)')\n",
    "axs[0].hist(r,bins='fd', color='k', alpha=.5, label='r')\n",
    "axs[0].legend(loc='upper right')\n",
    "axs[0].set_xlabel('Value')\n",
    "axs[0].set_ylabel('Frequency')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Histograms of r and F(r)')\n",
    "\n",
    "# scatter plot to visualize the transform\n",
    "axs[1].plot(r,fish_r,'ko',markersize=10,alpha=.5,markerfacecolor=(.8,.8,.8))\n",
    "axs[1].set_xlabel('r')\n",
    "axs[1].set_ylabel('F(r)')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Scatter plot of r vs. F(r)')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_fisherReal.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "jm8f6-zaS8Vs"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "2304rpsgK1VV"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Figure 12.11: Subgroups paradox"
   ],
   "metadata": {
    "id": "RKOSUH0LK1YG"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# initializations\n",
    "n = 20 # sample points per group\n",
    "offsets = [2, 3.5, 5] # mean offsets\n",
    "\n",
    "allx = np.array([])\n",
    "ally = np.array([])\n",
    "\n",
    "s = 'so^' # marker shapes\n",
    "\n",
    "# generate and plot data\n",
    "plt.figure(figsize=(5,6))\n",
    "for datai in range(3):\n",
    "\n",
    "  # generate data\n",
    "  x = np.linspace(offsets[datai]-1,offsets[datai]+1,n)\n",
    "  y = -x/3 + np.mean(x) + np.random.randn(n)/3\n",
    "\n",
    "  # subgroup correlation\n",
    "  r,p = stats.pearsonr(x,y)\n",
    "\n",
    "  # plot\n",
    "  plt.plot(x,y,s[datai],color=(datai/3,datai/3,datai/3),\n",
    "           markersize=14,alpha=.7,label=f'r={r:.2f}, p={p:.2f}')\n",
    "\n",
    "  # gather the data into one array\n",
    "  allx = np.append(allx,x)\n",
    "  ally = np.append(ally,y)\n",
    "\n",
    "\n",
    "\n",
    "# % now correlate the groups\n",
    "[r,p] = stats.pearsonr(allx,ally)\n",
    "plt.title(f'Total r={r:.2f}, p={p:.3f}',loc='center')\n",
    "\n",
    "plt.xticks([])\n",
    "plt.yticks([])\n",
    "plt.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_simpsons.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "RVjK-ozIMhrk"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "23F34ndkbhK0"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Cosine similarity"
   ],
   "metadata": {
    "id": "ggauNuUObhHr"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# variables\n",
    "x = np.array([ 1,2,3,4 ])\n",
    "y = np.array([ 1,2,3,4 ])\n",
    "z = np.array([ 101,102,103,104])\n",
    "\n",
    "# correlations\n",
    "r_xy = np.corrcoef(x,y)[0,1]\n",
    "r_xz = np.corrcoef(x,z)[0,1]\n",
    "\n",
    "# cosine similarities\n",
    "c_xy = 1-spatial.distance.cosine(x,y)\n",
    "c_xz = 1-spatial.distance.cosine(x,z)\n",
    "\n",
    "\n",
    "# print out the results\n",
    "print(f'r_xy: {r_xy:.3f}')\n",
    "print(f'c_xy: {c_xy:.3f}')\n",
    "print('')\n",
    "print(f'r_xz: {r_xz:.3f}')\n",
    "print(f'c_xz: {c_xz:.3f}')"
   ],
   "metadata": {
    "id": "-UwQ9baBbhFG"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "KPExLqIybQ57"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 1"
   ],
   "metadata": {
    "id": "FGXH90oCbQ3P"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# two random correlated variables\n",
    "v = np.random.randn(10)\n",
    "w = v + np.random.randn(len(v))\n",
    "\n",
    "### correlation using mean-centered dot products\n",
    "# mean-center\n",
    "vm = v-np.mean(v)\n",
    "wm = w-np.mean(w)\n",
    "\n",
    "# dot products\n",
    "r_me = np.dot(vm,wm) / np.sqrt(np.dot(vm,vm)*np.dot(wm,wm))\n",
    "\n",
    "\n",
    "### correlation using numpy\n",
    "r_np = np.corrcoef(v,w)[0,1]\n",
    "\n",
    "\n",
    "# print results\n",
    "print(f'r from np.corr: {r_np:.3f}')\n",
    "print(f'r from np.dot : {r_me:.3f}')"
   ],
   "metadata": {
    "id": "7lgkn1v-bSvz"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "KbetrRqTbSyh"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 2"
   ],
   "metadata": {
    "id": "W23OSbsucWy1"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create random data\n",
    "N = 43\n",
    "r = .4\n",
    "x = np.random.randn(N)\n",
    "y = np.random.randn(N)\n",
    "y = x*r + y*np.sqrt(1-r**2)\n",
    "\n",
    "# r,p from numpy\n",
    "r_np = np.corrcoef(x,y)[0,1]\n",
    "t = r_np*np.sqrt(N-2) / (1-r_np**2)\n",
    "p_np = stats.t.sf(np.abs(t),N-2) * 2  # times 2 for a two-sided test\n",
    "\n",
    "# r,p from scipy\n",
    "r_sp,p_sp = stats.pearsonr(x,y)\n",
    "\n",
    "# print correlation values\n",
    "print(f'r (p) from numpy: {r_np:.4f} ({p_np:.4f})')\n",
    "print(f'r (p) from scipy: {r_sp:.4f} ({p_sp:.4f})')"
   ],
   "metadata": {
    "id": "tIS-lV06f_l2"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# now in a loop\n",
    "\n",
    "# I wrote a function to output the two p-values,\n",
    "# which makes the experiment code below easier to read.\n",
    "def getpvals(x,y):\n",
    "\n",
    "  # r,p from numpy\n",
    "  r_np = np.corrcoef(x,y)[0,1]\n",
    "  t = r_np*np.sqrt(len(x)-2) / (1-r_np**2)\n",
    "  p_np = stats.t.sf(np.abs(t),len(x)-2) * 2  # times 2 for a two-sided test\n",
    "\n",
    "  # r,p from scipy\n",
    "  r_sp,p_sp = stats.pearsonr(x,y)\n",
    "\n",
    "  return p_np,p_sp\n",
    "\n",
    "\n",
    "# range of correlation values\n",
    "rvals = np.linspace(0,.99,40)\n",
    "\n",
    "# results matrix\n",
    "pvalues = np.zeros((len(rvals),2))\n",
    "\n",
    "\n",
    "## run the experiment\n",
    "for ri in range(len(rvals)):\n",
    "\n",
    "  # create the data\n",
    "  x = np.random.randn(44)\n",
    "  y = np.random.randn(44)\n",
    "  y = x*rvals[ri] + y*np.sqrt(1-rvals[ri]**2)\n",
    "\n",
    "  # get the two p-values\n",
    "  pvalues[ri,:] = getpvals(x,y)\n",
    "\n",
    "\n",
    "\n",
    "# plot\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.plot(rvals,np.log(pvalues[:,1]),'ko',markersize=12,markerfacecolor=(.9,.9,.9),label='scipy')\n",
    "plt.plot(rvals,np.log(pvalues[:,0]),'ks',markersize=12,markerfacecolor=(.7,.7,.7),label='numpy',zorder=-1)\n",
    "plt.legend()\n",
    "plt.xlabel('Correlation value (r)')\n",
    "plt.ylabel('log(p)')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_ex2.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "RPB96tMXcVRX"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "OaWjIAXncVOY"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 3"
   ],
   "metadata": {
    "id": "h1tKSPsubS1k"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# matrix of p-values\n",
    "\n",
    "N = 10000 # observations\n",
    "M = 15 # features\n",
    "\n",
    "# data matrix\n",
    "X = np.random.randn(N,M)\n",
    "\n",
    "# correlation matrix\n",
    "R = np.corrcoef(X.T)\n",
    "\n",
    "# confirm that it's the right shape\n",
    "print(f'Correlation matrix shape: {R.shape}')"
   ],
   "metadata": {
    "id": "bWE0DlYz-DMX"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# compute the t-values\n",
    "Tnum = R*np.sqrt(N-2)\n",
    "Tden = 1-R**2 + np.finfo(float).eps # adding a tiny number to avoid n/0\n",
    "\n",
    "T = Tnum / Tden\n",
    "\n",
    "# compute the p-values\n",
    "P = stats.t.sf(T,N-2)\n",
    "\n",
    "# visualize all matrices\n",
    "fig,axs = plt.subplots(1,3,figsize=(10,5))\n",
    "\n",
    "cax = axs[0].imshow(R,vmin=-.1,vmax=.1,cmap='gray')\n",
    "axs[0].set_title(r'$\\bf{A}$)  R matrix')\n",
    "c = fig.colorbar(cax,fraction=.046,pad=.04); c.ax.tick_params(labelsize=10)\n",
    "\n",
    "cax = axs[1].imshow(T,vmin=-2,vmax=2,cmap='gray')\n",
    "axs[1].set_title(r'$\\bf{B}$)  T matrix')\n",
    "c = fig.colorbar(cax,fraction=.046,pad=.04); c.ax.tick_params(labelsize=10)\n",
    "\n",
    "cax = axs[2].imshow(P,vmin=0,vmax=.05,cmap='gray')\n",
    "axs[2].set_title(r'$\\bf{C}$)  P matrix')\n",
    "c = fig.colorbar(cax,fraction=.046,pad=.04); c.ax.tick_params(labelsize=10)\n",
    "\n",
    "# properties common to all axes\n",
    "for a in axs:\n",
    "  a.set(xticks=[],yticks=[],xlabel='Features',ylabel='Features')\n",
    "  a.spines[['right','top']].set_visible(True)\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_ex3.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "tlLjNkq5-CAY"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "toH-Wa8g-B8p"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 4"
   ],
   "metadata": {
    "id": "3LqK0Fsd-B5t"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create the sigma matrix\n",
    "Sigma = np.std(X,ddof=1,axis=0)\n",
    "Sigma = np.diag(Sigma)\n",
    "\n",
    "# compute C from R\n",
    "C_me = Sigma@R@Sigma # from formula\n",
    "C_np = np.cov(X.T,ddof=1) # from numpy\n",
    "\n",
    "# check for equality (these should all be True)\n",
    "(C_np-C_me)"
   ],
   "metadata": {
    "id": "hdFig4LGIUrB"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Now compute R from C\n",
    "invSigma = 1/np.std(X,ddof=1,axis=0)\n",
    "invSigma = np.diag(Sigma)\n",
    "\n",
    "R_me = invSigma@C_np@invSigma\n",
    "\n",
    "\n",
    "# check for equality (these should all be True)\n",
    "(R-R_me) < 1e-16"
   ],
   "metadata": {
    "id": "9ildw1DgIUtd"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "0rUOGR_T-B2u"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 5"
   ],
   "metadata": {
    "id": "nrM1GNjRf_jG"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# simulation parameters\n",
    "N = 1000  # observations\n",
    "M =   20  # features\n",
    "\n",
    "# number of repetitions\n",
    "numReps = 30\n",
    "\n",
    "\n",
    "# initialize data lists\n",
    "alldata = np.zeros((M,N))\n",
    "corrmats = np.zeros((M,M,numReps+1))\n",
    "\n",
    "# \"pure\" data\n",
    "covars = np.linspace(-1,1,M)[:,None]\n",
    "dataOG = np.random.randn(N) * covars\n",
    "\n",
    "\n",
    "# random noise in each repetition\n",
    "for idx in range(numReps):\n",
    "\n",
    "  # this run's data\n",
    "  thisdata = dataOG + np.random.randn(M,N)*5\n",
    "\n",
    "  # its correlation matrix\n",
    "  corrmats[:,:,idx] = np.corrcoef(thisdata)\n",
    "\n",
    "  # sum the data\n",
    "  alldata += thisdata\n",
    "\n",
    "\n",
    "\n",
    "# correlation of data average\n",
    "corrmats[:,:,-1] = np.corrcoef(alldata)\n",
    "\n",
    "\n",
    "### plotting\n",
    "fig,axs = plt.subplots(1,3,figsize=(12,5))\n",
    "\n",
    "axs[0].imshow(covars@covars.T,vmin=-.3,vmax=.3,cmap='gray')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Ground truth')\n",
    "\n",
    "axs[1].imshow(np.mean(corrmats[:,:,:-1],axis=2),vmin=-.3,vmax=.3,cmap='gray')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Ave. of correlations')\n",
    "\n",
    "cax = axs[2].imshow(corrmats[:,:,-1],vmin=-.3,vmax=.3,cmap='gray')\n",
    "axs[2].set_title(r'$\\bf{C}$)  Correlation of ave.')\n",
    "\n",
    "cbar_ax = fig.add_axes([.91,.22,.015,.55])\n",
    "cbar = plt.colorbar(cax,cax=cbar_ax)\n",
    "\n",
    "\n",
    "# properties common to all axes\n",
    "for a in axs:\n",
    "  a.set(xticks=[],yticks=[],xlabel='Features',ylabel='Features')\n",
    "  a.spines[['right','top']].set_visible(True)\n",
    "\n",
    "# plt.tight_layout()\n",
    "plt.savefig('cor_ex5.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "Przdhj9bvmAT"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "_d-rCnvvfXGs"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 6"
   ],
   "metadata": {
    "id": "4duacQH146c_"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# re-initialize data lists\n",
    "alldata = np.zeros((M,N))\n",
    "\n",
    "\n",
    "### run the experiment\n",
    "for idx in range(numReps):\n",
    "\n",
    "  # this run's data (only 'covars' is constant across repetitions)\n",
    "  thisdata = np.random.randn(N)*covars + np.random.randn(M,N)\n",
    "\n",
    "  # its correlation matrix\n",
    "  corrmats[:,:,idx] = np.corrcoef(thisdata)\n",
    "\n",
    "  # store the data\n",
    "  alldata += thisdata\n",
    "\n",
    "\n",
    "\n",
    "# correlation of data average\n",
    "corrmats[:,:,-1] = np.corrcoef(alldata)\n",
    "\n",
    "\n",
    "### plotting\n",
    "fig,axs = plt.subplots(1,3,figsize=(12,5))\n",
    "\n",
    "axs[0].imshow(covars@covars.T,vmin=-.3,vmax=.3,cmap='gray')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Ground truth')\n",
    "\n",
    "axs[1].imshow(np.mean(corrmats[:,:,:-1],axis=2),vmin=-.3,vmax=.3,cmap='gray')\n",
    "axs[1].set_title(r'$\\bf{B}$)  Ave. of correlations')\n",
    "\n",
    "cax = axs[2].imshow(corrmats[:,:,-1],vmin=-.3,vmax=.3,cmap='gray')\n",
    "axs[2].set_title(r'$\\bf{C}$)  Correlation of ave.')\n",
    "\n",
    "cbar_ax = fig.add_axes([.91,.22,.015,.55])\n",
    "cbar = plt.colorbar(cax,cax=cbar_ax)\n",
    "\n",
    "\n",
    "# properties common to all axes\n",
    "for a in axs:\n",
    "  a.set(xticks=[],yticks=[],xlabel='Features',ylabel='Features')\n",
    "  a.spines[['right','top']].set_visible(True)\n",
    "\n",
    "# plt.tight_layout()\n",
    "plt.savefig('cor_ex6.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "QugiU8PKf_gP"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "hbReUZ-D-Bzw"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 7"
   ],
   "metadata": {
    "id": "zNANilrvbS90"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "samplesize = 30\n",
    "nSamples = 23\n",
    "\n",
    "corrs = np.zeros(nSamples)\n",
    "tres  = np.zeros(2)\n",
    "\n",
    "# loop over experiments\n",
    "for ni in range(nSamples):\n",
    "\n",
    "  # create the data\n",
    "  x = np.random.randn(samplesize)\n",
    "  y = np.random.randn(samplesize)\n",
    "  y = x*.1 + y*np.sqrt(1-.1**2)\n",
    "\n",
    "  # correlation\n",
    "  corrs[ni] = np.corrcoef(x,y)[0,1]\n",
    "\n",
    "# now for a t-test on r values\n",
    "tres[0] = stats.ttest_1samp(corrs,0).statistic\n",
    "tres[1] = stats.ttest_1samp(np.arctanh(corrs),0).statistic\n",
    "\n",
    "# critical t-value\n",
    "tCrit = stats.t.isf(.05/2,samplesize-2)\n",
    "\n",
    "\n",
    "print(f't-value from \"raw\" coefficients   : {tres[0]:.4f}')\n",
    "print(f't-value from Fisher-z coefficients: {tres[1]:.4f}')\n",
    "print(f'Critical t-value for p<.05        : {tCrit:.4f}')"
   ],
   "metadata": {
    "id": "GkfxnjI3D5cA"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# now for a wider range of r values\n",
    "rs = np.linspace(.01,.5,21)\n",
    "samplesize = 30\n",
    "nSamples = 23\n",
    "\n",
    "corrs = np.zeros((len(rs),nSamples))\n",
    "tres  = np.zeros((len(rs),2))\n",
    "\n",
    "# critical t (doesn't depend on the population r or sample size)\n",
    "tCrit = stats.t.isf(.05/2,samplesize-2)\n",
    "\n",
    "\n",
    "# run the experiment!\n",
    "for ri,r in enumerate(rs):\n",
    "\n",
    "  # loop over experiments\n",
    "  for ni in range(nSamples):\n",
    "\n",
    "    # create the data\n",
    "    x = np.random.randn(samplesize)\n",
    "    y = np.random.randn(samplesize)\n",
    "    y = x*r + y*np.sqrt(1-r**2)\n",
    "\n",
    "    # correlation\n",
    "    corrs[ri,ni] = np.corrcoef(x,y)[0,1]\n",
    "\n",
    "  # now for a t-test on r values\n",
    "  tres[ri,0] = stats.ttest_1samp(corrs[ri,:],0).statistic\n",
    "  tres[ri,1] = stats.ttest_1samp(np.arctanh(corrs[ri,:]),0).statistic\n",
    "\n",
    "\n",
    "\n",
    "## plot\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4))\n",
    "axs[0].plot(rs,np.mean(corrs,axis=1),'ks',markersize=10,markerfacecolor=(.6,.6,.6))\n",
    "axs[0].plot([rs[0],rs[-1]],[rs[0],rs[-1]],'--',color=(.8,.8,.8),zorder=-3)\n",
    "axs[0].set(xlabel='Theoretical correlation',ylabel='Empirical correlation')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Theoretical vs. empirical r')\n",
    "\n",
    "axs[1].plot(rs,tres[:,0],'ks',markersize=10,markerfacecolor=(.6,.6,.6),label='\"Raw\"')\n",
    "axs[1].plot(rs,tres[:,1],'ko',markersize=10,markerfacecolor=(.3,.3,.3),label='Fisher-z')\n",
    "axs[1].plot(rs,np.full(len(rs),tCrit),'--',color=(.8,.8,.8),zorder=-3,label='Critical t')\n",
    "axs[1].legend()\n",
    "axs[1].set(xlabel='Theoretical correlation',ylabel='T values')\n",
    "axs[1].set_title(r\"$\\bf{B}$)  T-values ($r's\\neq 0$)\")\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_ex7.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "B0x7wpf_f_ox"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "CC2CrtLZcW1P"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 8"
   ],
   "metadata": {
    "id": "wFH0I4B8cVLj"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# generate some correlated random data\n",
    "x = np.random.randn(40)\n",
    "y = x + np.random.randn(len(x))\n",
    "\n",
    "# manual cosine similarity\n",
    "cs_num = sum(x*y)\n",
    "cs_den = np.sqrt(sum(x*x)) * np.sqrt(sum(y*y))\n",
    "cs_me = cs_num / cs_den\n",
    "\n",
    "# using the  distance function in the scipy.spatial library\n",
    "# Note: using this function is confusing, because it computes *distance* although we want *similarity*.\n",
    "# Fortunately, the two are simple inverses, so one is 1- the other.\n",
    "cs_sp = 1-spatial.distance.cosine(x,y)\n",
    "\n",
    "\n",
    "print(f'Manual result: {cs_me:.3f}')\n",
    "print(f'Scipy.spatial: {cs_sp:.3f}')"
   ],
   "metadata": {
    "id": "OdX0bsZAQkaq"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "DZeQpfsURHae"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 9"
   ],
   "metadata": {
    "id": "ljIxhxvNf_Xk"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "#  range of requested correlation coefficients\n",
    "rs = np.linspace(-1,1,100)\n",
    "\n",
    "# sample size\n",
    "N = 500\n",
    "\n",
    "\n",
    "# initialize output matrix\n",
    "corrs = np.zeros((len(rs),2))\n",
    "\n",
    "\n",
    "# loop over a range of r values\n",
    "for ri in range(len(rs)):\n",
    "\n",
    "  # generate data\n",
    "  x = np.random.randn(N)\n",
    "  y = x*rs[ri] + np.random.randn(N)*np.sqrt(1-rs[ri]**2)\n",
    "\n",
    "  # mean de-centering\n",
    "  x = x-10\n",
    "\n",
    "  # compute correlation\n",
    "  corrs[ri,0] = np.corrcoef(x,y)[0,1]\n",
    "\n",
    "  # compute cosine similarity\n",
    "  corrs[ri,1] = 1-spatial.distance.cosine(x,y)\n",
    "\n",
    "\n",
    "\n",
    "## visualize the results\n",
    "_,axs = plt.subplots(1,2,figsize=(10,4.5))\n",
    "\n",
    "axs[0].plot(rs,corrs[:,0],'ks',markersize=10,markerfacecolor=(.5,.5,.5),alpha=.5,label='Correlation')\n",
    "axs[0].plot(rs,corrs[:,1],'ko',markersize=10,markerfacecolor=(.9,.9,.9),alpha=.5,label='Cosine sim.')\n",
    "axs[0].legend()\n",
    "axs[0].set(xlabel='Requested correlation',ylabel=r'Empirical $r$ or $S_C$')\n",
    "axs[0].set_title(r'$\\bf{A}$)  Correlation and cosine sim.')\n",
    "\n",
    "axs[1].plot(corrs[:,0],corrs[:,1],'ks',markersize=10,markerfacecolor=(.2,.2,.2),alpha=.5)\n",
    "axs[1].axhline(y=0,color='gray',linestyle='--')\n",
    "axs[1].axvline(x=0,color='gray',linestyle='--')\n",
    "axs[1].set(xlabel='Correlation',ylabel='Cosine similarity')\n",
    "axs[1].set_title(rf'$\\bf{{B}}$)  r={np.corrcoef(corrs.T)[1,0]:.2f}')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_ex9.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "cPtol2DcOk0d"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "9pHlVQOXf_U1"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 10"
   ],
   "metadata": {
    "id": "2wK-0xHGbQ0P"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# import the data\n",
    "url = \"https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data\"\n",
    "column_names = ['MPG','Cylinders','Displacement','Horsepower','Weight','Acceleration','Model Year','Origin','Car Name']\n",
    "\n",
    "data = pd.read_csv(url,delim_whitespace=True,names=column_names, na_values=\"?\")\n",
    "data"
   ],
   "metadata": {
    "id": "XRbMROryeYNd"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# examine distributions\n",
    "\n",
    "# include only numerical variables\n",
    "data_numerical = data.drop(columns=['Car Name','Origin'])\n",
    "\n",
    "# draw histograms with seaborn\n",
    "fig,axs = plt.subplots(2,4,figsize=(12,6))\n",
    "for a,column in zip(axs.flatten(),data_numerical.columns):\n",
    "  sns.histplot(data=data_numerical, x=column, ax=a, color=(.7,.7,.7))\n",
    "\n",
    "# only 7 columns, so switch off the empty 8th :P\n",
    "axs[-1,-1].axis('off')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_ex10b.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "SANu3BsYeYP6"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Create a correlation matrix\n",
    "R = data_numerical.corr(method='spearman')\n",
    "\n",
    "plt.figure(figsize=(8,6))\n",
    "sns.heatmap(R, annot=True, cmap='coolwarm',vmin=-1,\n",
    "            xticklabels=R.columns,yticklabels=R.columns)\n",
    "\n",
    "plt.title('Correlation matrix (Spearman)',loc='center',weight='bold')\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_ex10c.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "aBoJS_pJmnww"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "LqYEQR1FZGYr"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 11"
   ],
   "metadata": {
    "id": "O1xbD8nSZGQW"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Calculate R and P matrices from Spearman correlation\n",
    "R,P = stats.spearmanr(data_numerical,nan_policy='omit')\n",
    "\n",
    "# store as dataframes for seaborn plotting (note: \"df\" here means \"dataframe\")\n",
    "R_df = pd.DataFrame(R, columns=data_numerical.columns, index=data_numerical.columns)\n",
    "P_df = pd.DataFrame(P, columns=data_numerical.columns, index=data_numerical.columns)\n",
    "\n",
    "# Bonferroni correction [ formula is (M*(M-1))/2 ]\n",
    "num_comparisons = (data_numerical.shape[1]*(data_numerical.shape[1]-1)) / 2\n",
    "bonferroni_thresh = .05 / num_comparisons\n",
    "significant = P_df < bonferroni_thresh\n",
    "\n",
    "\n",
    "# Create a matrix of annotations\n",
    "annot_array = R_df.astype(str).values\n",
    "\n",
    "# loop through all elements of the matrix and create a string to display\n",
    "for i in range(R_df.shape[0]):\n",
    "  for j in range(R_df.shape[1]):\n",
    "\n",
    "    # the string depends on the significance\n",
    "    if not significant.iloc[i,j]:\n",
    "      # if non-significant, just the correlation coefficient\n",
    "      annot_array[i,j] = f'{R_df.iloc[i, j]:.2f}'\n",
    "    else:\n",
    "      # if significant, add an asterisk to the coefficient\n",
    "      annot_array[i,j] = f'{R_df.iloc[i, j]:.2f}*'\n",
    "\n",
    "    # don't need to report the diagonals (trivially=1)\n",
    "    if i==j:\n",
    "      annot_array[i,j] = ''\n",
    "\n",
    "\n",
    "## now show the image\n",
    "plt.figure(figsize=(8,6))\n",
    "sns.heatmap(R_df,annot=annot_array,fmt='s',cmap='coolwarm',vmin=-1,\n",
    "            xticklabels=R_df.columns,yticklabels=R_df.columns)\n",
    "\n",
    "plt.title('Correlation matrix (*p<.05 corrected)',loc='center',weight='bold')\n",
    "plt.tight_layout()\n",
    "plt.savefig('cor_ex11.png')\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "Pd2_gg0feYY-"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "id": "JutlK7jseYs6"
   },
   "execution_count": null,
   "outputs": []
  }
 ]
}